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We present a method based on the Path Integral Monte Carlo formalism for the 
calculation of ground-state time correlation functions in quantum systems. The key 
point of the method is the consideration of time as a complex variable whose phase 5 
acts as an adjustable parameter. By using high-order approximations for the quan¬ 
tum propagator, it is possible to obtain Monte Carlo data all the way from purely 
imaginary time to 5 values near the limit of real time. As a consequence, it is possible 
to infer accurately the spectral functions using simple inversion algorithms. We test 
this approach in the calculation of the dynamic structure function S{q,uj) of two 
one-dimensional model systems, harmonic and quartic oscillators, for which S{q,u) 
can be exactly calculated. We notice a clear improvement in the calculation of the 
dynamic response with respect to the common approach based on the inverse Laplace 
transform of the imaginary-time correlation function. 
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I. INTRODUCTION 


In the last decades, quantum Monte Carlo (QMC) methods have been extensively used 
in the held of quantum many-body physics. Many of these numerical techniques rely on 
stochastic propagation in imaginary time and can provide extremely accurate results for the 
thermodynamic and static properties of many-body systems, even in those where quantum 
correlations make unavoidable the use of non-perturbative approaches.-*- The main draw¬ 
back of QMC methods is the difficulty arising in the calculation of spectral functions. These 
functions, which are particularly relevant for the study of the dynamical properties of quan¬ 
tum many-body systems (e.g. the excitation spectrum or the transport coefficients), can 
be obtained as Fourier transforms of real-time correlation functions. A QMC calculation 
of these quantities, however, is particularly inefficient since the rapidly oscillating exponen¬ 
tials appearing in the dehnition of real-time propagators make the statistical errors grow 
exponentially with time. Many approximation schemes have been developed and used to 
investigate the dynamic properties of quantum many-body systems. For instance, centroid^ 
or ring-polymer molecular dynamics^ has been successfully applied to the study of quan¬ 
tum many-body systems in the semi-classical regime. In the limit of zero temperature, an 
alternative approach is to use correlated perturbation theory^ relying on the ground-state 
properties of the system obtained with QMC calculations.-*— 

Nevertheless, the mainstream approaches to the calculation of spectral functions from 
QMC simulations consist in attempting a numerical inversion of a Laplace transform. This 
integral transform relates the desired spectral functions to the correlation functions in imag¬ 
inary time, easily attainable with QMC methods. However, the inverse Laplace transform 
of noisy data is an ill-posed problem. This means that, given a particular set of data for the 
imaginary-time correlation function, it is hardly possible to recover a unique, well-dehned 
solution to the problem. Sophisticated regularization techniques can then be used to repro¬ 
duce a reasonable estimate of the spectral function.— In the last decades, several algorithms 
to deal with the inverse Laplace transform of noisy data have been proposed,—*— but these 
methods can only be reliably applied to the analysis of the low-energy dynamic properties 
of quantum systems, since the Laplace kernel tends to suppress high-energy contributions. 
In order to overcome these limitations and to get more accurate results of spectral functions 
from QMC data, it is necessary to develop new estimators for the quantum time correlation 
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functions. — 

In this work, we propose to infer the dynamic structure function of a quantum system at 
zero temperature from a QMC estimation of the corresponding correlation function in com¬ 
plex time. Similar approaches have been already used for studying the dynamic properties 
of quantum systems at finite temperature T. In this case, the term (with (3 = 1/T) 
appearing in the definition of the thermal averages can be considered as an evolution oper¬ 
ator in imaginary time. Thus, the real-time correlation function can be rewritten in terms 
of a correlation function in complex time,— which can be calculated using path-integral 
formalisms^ and estimated in QMC calculations.— Even though this estimation is reliable 
only for times t < h(3, the spectral functions obtained within this approach exhibit a sig- 
nihcant improvement over the results derived from analytic continuation of imaginary-time 
correlation functions.—^— 

Our goal is to extend this formalism to the calculation of ground-state time correlation 
functions, even considering that at zero temperature the notion of complex time has not 
a precise physical meaning. This strategy allows us to introduce an adjustable parameter, 
namely the phase S of the complex time tc = which makes possible to calculate the 

correlation function in an intermediate regime between the commonly used imaginary time 
{S = 7^/2) and real time (5 = 0). 

More precisely, we sample paths connecting two configurations distributed according to 
the ground-state wave function of the quantum system and calculate, over these paths, the 
propagator at the complex time tc- Changing the phase 5, we can find an optimal value for 
which the correlation functions estimated with QMC are affected by moderate statistical 
errors and, at the same time, present a relevant amount of information on the real dynamics 
of the quantum system. This approach makes it possible to infer the spectral functions using 
rather simple inversion techniques since the ill-posed character of the inversion procedure 
is appreciably reduced. In this way, more accurate and more stable results than the usual 
ones, based on the inverse Laplace transform of imaginary-time data, can be obtained. 

Similarly to what happens in the case at finite temperature, the QMC estimation of the 
ground state correlation function in complex time is reliable only up to a certain value of |tc| 
depending on 5, above which the statistical error becomes too large and makes the numerical 
results meaningless. It is therefore crucial to develop strategies that make the range of 
accessible times as large as possible. In this work, we propose to tackle this problem using 
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high-order approximations for the quantum propagator.— In particular, we show that the 
propagator derived by Zillich et alM- is particularly suitable for the complex-time evolution. 

The rest of the paper is organized as follows. In Sec. II, we discuss the QMC method 
that we have devised for the calculation of complex-time correlation functions and, more 
briefly, the inversion method that we have used to obtain the dynamic structure function. In 
Sec. Ill, the results obtained with this method in one-dimensional problems are shown and 
compared with the standard approach relying only on imaginary-time correlation functions. 
Finally, the summary and main conclusions are reported in Sec. IV. 

II. METHOD 

A. Calculation of the complex-time correlation function 

The main objective of our work is the calculation of ground-state time correlation func¬ 
tions of a quantum system. At zero temperature, a general time correlation function is 
defined as 

CABitc) = , ( 1 ) 

where A and B are time-independent quantum mechanical operators in the Schrodinger 
picture corresponding to measurable observables, H is the Hamiltonian, and |\l/o) is the 
ground state. For the sake of simplicity, in the following we consider correlations among 
operators which are diagonal in coordinate space and use one-dimensional notation (the 
generalization to multi-dimensional space is straightforward). 

The main idea of this work is to calculate CAB{tc), defined in Eq. [H where tc has been 
analytically extended to the complex plane. We indicate with > 0 and —6 the modulus 
and the phase of the complex time, tc = respectively. In order to elaborate a form 

for the estimator of CAB{tc) implementable in computer simulations, we rewrite Eq. [T] in 
the coordinate space, 

CAB{tc) = J dxodxMe'^‘=^°{^o\xM){xM\Ae~'^‘^^B\xo){xo\'^o) = 

= -^ J dxodxM%{xM)A{xM)G{xo,XM;tc)B{xo)^o{xo) , (2) 

where G{xo,XM',tc) = {xm\^~^^''^\xo) is the propagator from position xq to position xm in 
complex time tc, and A/” is a normalization constant. In the general case of complex time 
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tc, the propagator G{xo,XM',tc) is a complex function that becomes real and positive only 
when tc is a purely imaginary time. Thus, the function G{xQ,XM]tc) cannot be used as a 
probability distribution function for the sampling of coordinates in any QMC algorithm (as 
it is normally done, for instance, in the PIMC method). Therefore, what we do is to sample 
hrst the positions xq and xm according to a probability distribution constructed from an 
accurate approximation to the ground-state wave function. This sampling can be performed 
using any conventional QMC technique at zero temperature. In this work, we use the Path 
Integral Ground State (PIGS) method.-*^ Then, having sampled the positions Xq and xm, 
we calculate GAsitc) estimating the quantity A{xM)G{xo,XM]tc)B{xo). 

In order to carry on this procedure one needs to know the exact form of the Green’s 
function G(xo,xm; 4) for any value but this is in general unknown. However, what 
is possible is to construct accurate approximations to the propagator in the limit of small 
tm = \tc\- Then, to estimate GABitc) for larger values of tm we use the path-integral formalism 
to rewrite G{xo,XM]tc) as a convolution of M propagators of a shorter time Ec = tc/M, 

„ M 

G{xo,XM]tc) = / dxi -. .dxM-iY[G{xk,Xk-i]ec) . (3) 

k=i 

Within this approach, it becomes necessary to sample all the conhgurations {xi,X2, ■ ■ ■, xm-i}, 
i.e., to build paths from the position xq to the position xm- However, the choice of the 
probability distribution Ppath(3^0) • • • ,xm) for these paths is not trivial and depends on 
the system studied. Generally, we notice that using imaginary-time propagator to this end 
is not a good choice, because in this case the sampled paths would remain close to the 
minimum energy path and the estimator would not be able to capture all the contributions 
to Gab coming from the excited states. As a simple and flexible enough option, it is possible 
to choose Ppath as the product of M free propagators of imaginary-time step r^. 


M 


Ppath (^0: ; • 

. . . , Xm) GfYeei^Xf^^ T^) , 

k=l 

( 4 ) 

G^free(^/i;5 ^k—1^ Ts) 

-(47rAT,)“-'"exp(' ' 

( 5 ) 


In Eq. m is the number of particles, d is the dimensionality of the system, and A = 
h?/{2m). This choice indeed allows to construct the paths by means of simple sampling 
techniques which do not require a large computational effort, like for instance the staging 
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algorithm.— In the case of quantum systems interacting with a smooth potential, we 
notice that it is possible to obtain good results for CABiic) using Ppath iu Eq. |H provided 
that the parameter r* is properly chosen. Indeed, we see that the variance of the estimator 
for CAsitc) is reduced when the free propagator in the imaginary time r* is similar to the 
modulus of the kinetic propagator in the complex time Sc- 

Since the purpose of this work is to test our QMC approach in two model systems inter¬ 
acting with smooth potentials (the quantum harmonic and quartic oscillators), we decide 
to use this choice of Ppath with ~ to perform the sampling of the paths 

{a:i,X 2 ,... Nevertheless, this may not be the best choice in general, and one may 

have to use more sophisticated and more computationally demanding algorithms for the 
sampling of the paths. 

Once the probability distribution Ppath is chosen, the expression of the ground state 
complex time correlation function becomes 

CAsitc) = AA' [ dxQ... dxMA{xM) ^''\ b{xq) x 

J Ppathi^O? • • • 5 ^m) 

'^o{xM)Ppi,th{xo, Xi, ..., XM)d^o(xo) • (6) 

At this point, one has to choose an approximation scheme for G{xk,Xk-i;Sc) in order to 
derive an analytical expression that can be implemented in computer simulations. Increasing 
the number of convolution terms M, and thus decreasing the modulus of Ec, it is possible 
to systematically improve the quality of the approximation and to asymptotically recover 
the exact correlation function. Nevertheless, every propagator G{xk,Xk-i;£c) introduces an 
oscillating phase term in the integrand of Eq. [6], and thus the statistical noise of the estimator 
for GAsitc) increases notably when M becomes large. In order to obtain reliable results, 
it is fundamental to develop numerical strategies that keep the number M of convolution 
terms as low as possible. 

The simplest approximation to the propagator is the primitive approximation (PA), which 
relies on the factorization ~ where K and V are the kinetic and potential 

operators, respectively. In this scheme, the complex-time propagator can be written as 


G{Xki Xk—li E(^ — GpAiXki Xk—li 

(Xk - Xfc-i)^ 


= exp 


4A i£r 


exp 


.V(xk) + V(xk-i) 
2h 


(7) 


The PA approximation is easily implementable within our QMC procedure but requires a 
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large number M of convolution terms in Eq. El In order to improve the accuracy, it is impor¬ 
tant to use higher-order approximations to the complex-time propagator. In conventional 
PIMC simulations, a signihcant improvement in efficiency can be obtained using symplectic 
expansions of the time-evolution operator that incorporates double commutators between 
kinetic and potential operators.— >2^ For local potentials, these commutators lead to extra 
terms that are exponentials of the gradient of the potential squared times the third power 
of the time step. The inclusion of this contribution in the propagator largely improves the 
efficiency of the PIMC— and PIGS^ methods. In imaginary-time propagation, the contri¬ 
bution of the double commutator always appears in the argument of the exponential with 
a negative sign. However, in complex time this sign turns out to be positive for 5 < 60°, 
producing largely increasing amplitudes and thus unreliable results that make the use of this 
high-order scheme unpractical (see Appendix A). Therefore, it is very important to look for 
other expansions which can improve the PA but that do not include double-commutator 
terms. 

A high-order approximation for the complex-time propagator without double commuta¬ 
tor has been reported in Ref. In that work, the authors were able to improve the quality 
of the small-time propagator by introducing a linear combination, with some negative co¬ 
efficients, of different symplectic expansions on the same time. This expansion has some 
drawbacks when used in conventional PIMC simulations, since it gives rise to an approxi¬ 
mation for the imaginary-time propagator which is not positive dehnite. This feature does 
not represent a problem here, since in the calculation of CAsiic) the complex-time propa¬ 
gator is not used as the probability distribution of the Monte Carlo sampling but rather as 
the estimator. 

Once we have chosen the approximation for the complex-time propagator, the only thing 
that is still lacking in order to calculate CAB{tc) is the normalization constant J\f'. This can 
be computed imposing the autocorrelation function of the identity operator to be 1 for any 
value of tc- Therefore, if we dehne the complex quantity 

Kd/i;— ^^k—1^ ^c) 


Oa{Xo, . . . , Xm) — 


( 8 ) 


Ppath(^O) Xi, . . . , X ]\/[) 

where GA{xk,Xk-i',Ec) is the chosen approximation for the time propagator, the complex- 
time correlation function in Eq. E] can be written as 

{A{xm)Oa{xo, ..., xm)B{xo)) 


Cab(A) — 


{Oa{xo, .. .,Xm)) 


(9) 
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The bracket (...) indicates the averages over the conhgurations {xo,Xi,... ,xm} sampled 
following the scheme described above, i.e., with xq and xm sampled according to a reasonable 
approximation of the ground-state wave function, and {xi, X 2 , ■ ■ ■, xm-i} sampled according 
to the probability distribution Ppath{xo, xi,, xm)- 

Summarizing, the evaluation of CAsitc) © for a given complex time tc = consists 

of the following steps: 

1. To generate the Xq and xm configurations according to the probability distribution 
To(a;o)tko(a;M), by means of a QMC technique at zero temperature, like the PIGS 
algorithm. 

2. To choose M (number of points of the discrete path from xq to xm), so that the 

parameter Sm = is sufficiently small to recover the Sm ^ 0 limit. In practice, 

one selects the value of M that makes Em = tm/M < where the parameter 
depends on the accuracy of the approximated action. 

3. To generate the configurations {a;i,X 2 ,... ,xm-i}, he., the path from xq to xm, ac¬ 
cording to the probability distribution Ppath- 

4. To evaluate Oa{xo, ... ,xm) from Eq. [8] and accumulate the estimator of CAsitc) 
defined in Eq. [9l 

B. Inversion technique 

Once we have obtained the QMC data for the complex-time correlation function CAsitc), 
we need to recover the desired spectral function Sab{^) inverting the integral transform 



( 10 ) 


Considering that both the function CAsitc) and Sab{(x) are evaluated over a finite set of 


complex times {td} and frequencies {cujj, Eq. [10] is formally equivalent to a linear equation 


y = Ax , 


( 11 ) 


where the vector y represents the QMC data for the correlation function CABitc), the vector 
X the spectral function Sab{^) that we want to obtain, and A is a matrix defined from the 


kernel of the integral transform IfTUll which relates CAsitc) and Sab{^)- Notice that CAsitc) 
is a complex function: thus, its real and its imaginary part provide two different rows of the 
matrix A, both of them real. 

The best least-squares solution to Eq. [TT] is given by the pseudo-inverse matrix 

X = A^{AA^)-^y . (12) 

In well-posed problems, Eq. [12] directly provides useful solutions. If x has larger dimen¬ 
sionality than y, then the linear equation in 111 111 has an inhnite number of solutions, and 
flT^ provides the one which minimizes \x\^. Contrarily, if x has lower dimensionality than 
y, then no solution exists and Eq. [12] (using the Moore-Penrose pseudoinverse if AAA is not 
full-rank) provides the x vector which minimizes \y — Axp, i.e., a best ht to the y data is 
obtained. 

However, when the eigenvalues of the matrix AA^, which are all positive or zero, span 
a range of many orders of magnitude (in the numerical inversions performed in the present 
work, eigenvalues of covering the range 10°-10“^° are routinely found), the inversion 
problem becomes ill-posed, and the solution x to Eq. [12] is extremely sensitive to errors in 
the vector y. The ill-posed nature of the inversion process means that the statistical noise 
in the original data for C'As(tc), that is unavoidable in any QMC calculation, is uncontrol¬ 
lably magnified in the inversion process, resulting in a meaningless solution for the spectral 
function SabA)- 

In these situations, regularization techniques are useful to obtain meaningful solutions 
to the ill-posed problem.— The basic idea of these methods is to define a well-conditioned 
linear operator Ca which depends on a regularization parameter a > 0 that approaches the 
pseudo-inverse = A^{AA"^)~^ in the limit a —)■ 0. Then, the solution of the original 
problem can be obtained as x = hma_j.o Cay- 

In this work, we have chosen to use the Tikhonov regularization,— in which 

Ca = A^{AA^+ lA)-\ (13) 

where I is the identity matrix. Thanks to Tikhonov regularization, the solution x of the 
problem is much less sensitive to errors in the initial vector y. On the other hand, the 
regularization procedure introduces a bias in the estimation of x. The goal is however to 
keep the regularization parameter a as small as possible yo avoid introducing unwanted 
artifacts in the reconstructed solution. 
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In practice, the choice of the regularization parameter must avoid two different problems. 
If the regularization parameter a is too small, the solution is unstable and similar QMC data 
for the correlation function lead to different spectral functions. If a is too large, systematic 
effects start to appear in the solution. These effects can be controlled verifying that the 
correlation function obtained applying the direct integral transform (Eq. fTOl) to the given 
solution for the spectral function is in agreement with the starting QMC data for C^s(fc) 
(see Appendix B for additional information). Monte Carlo data of higher quality allow for 
smaller values of the regularization parameter and thus they are crucial for a satisfactory 
direct inversion. 

Focusing on the dynamic structure factor, the physical solution must verify Xi > 0 for 
every component of x since S{q,uj) > 0 . We introduce this requirement explicitly in the 
construction of the solution, making use of a square diagonal matrix Q = Diag(gi,..., q^), 
where each of the qi is to be understood as a factor (which we restrict to be either 0 or 1) 
that will multiply explicitly the component Xi of the vector solution x. The new solution, 
that can be written formally as 

X = Q A^{AQ A^)~^ y , (14) 

satishes by construction both Xj = 0 if g* = 0 and y = Ax, irrespective of Q. The reg¬ 
ularization procedure can be performed in this case by simply making the substitution 
A Q A^ A Q AA -I -1 a?. We use Eq. [Il]as a means of imposing the positiveness of S{q, cu). 
In order to do so, we set an iterative procedure, starting with Q = Diag(gi = 1,..., gAr = 1), 
using the regularized version of Eq. [TTl to obtain the vector solution x, and we set g* = 0 for 
all components Xj < 0 and form a new Q matrix which contains more zeroes in the diagonal 
than the previous one. Inserting the new Q back in Eq. [TTl a new solution is obtained. 
The procedure is repeated until no negative components are present, and we end up with a 
regularized, positive solution to the inversion problem. 

III. RESULTS 

The formalism developed in Sec. II has been applied to the calculation of the density- 
density correlation function in complex time, 

S{q,Q = , (15) 
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FIG. 1. (Color online) Real (top) and imaginary (bottom) parts of S{q,tc) for HP, with q = 1.5 
and 5 = as a function of tm = |tc|- The line stands for the exact result (fT6]l and the points 
to different approximations for the action. Triangles, primitive action; squares. Chin action^^ 
diamonds, Zillich action.— 


with the density-fluctuation operator pq = complex time tc = ■ The 

reliability of the method has been checked in two model problems which can be easily 
solved: a particle in a one-dimensional harmonic potential (HP), V(x) = x^/2, and a particle 
in a one-dimensional anharmonic potential (AP), V(x) = x^/4. We work in units where 
h = m = 1. The ground-state wave function Tq ffTHll is obtained using the PIGS algorithm 
with the high-order Chin action.— 

As commented in Sec. II, a relevant aspect that makes the calculation in complex time 
be more accurate is to use high-order actions in the evaluation of Eq. [TSl We need to work 
with as few number of beads M as possible to reach the maximum accessible time. In Fig. 
m we show results for the real and imaginary parts of S'(g, tc) for the HP as a function of tm. 
The results correspond to g = 1.5 and 5 = vr/Q. The line stands for the HP exact result,— 


S{q,tc) = exp 




1 ) 


(16) 


In the figure, we compare the exact function ffTHj) with our QMC results obtained with a 
single bead, M = 1, using different approximations for the actions employed in the evaluation 
of S{q,tc). As expected, the PA is only accurate at very short times. If we consider QMC 
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results for S{q,tc) with a relative error of 0.4%, we notice that these are in agreement with 
the exact result for tm ^ 0.3 and depart signihcantly of the exact result al larger time. 
Therefore, the PA is not a good choice because we would need a large number of beads to 
span the full time range. The results are signihcantly better if one uses high-order actions. 
In the hgure, we show estimations of the real and imaginary parts using the Chin action^ 
and a sixth-order expansion reported by Zillich et alM- Comparing numerical results of 
S{q,tc), with the same precision as before, we notice that the Chin action reproduces the 
exact results up to — 2. However, the Chin action is in general not appropriate because 
of the divergence terms derived from the double commutator (notice that for the HP this 
divergence is reduced because this contribution produces a renormalization of the oscillator 
frequency). The best result is obtained using the sixth-order approximation.— This action is 
able to account for the exact data up to tm — 3.5 and with the added beneht of not requiring 
double-commutator terms since it is based on extrapolations of PA actions with different 
time steps. Therefore, we have selected this action as the best option for this complex-time 
estimation. 

A second step in our methodology is the estimation of (see Sec. H) which determines 
the maximum time tm that can be covered with a single bead, with no signihcant bias coming 
from the small-time approximation of the action. This estimation is performed by studying 
the convergence of S{q,tc), with tm = \tc\ fixed, for small values of Em = tm/M. To perform 
this analysis, we have selected 5 = 7r/2 (imaginary time). Using a different value of 5, the 
statistical error of S{q,tc) tends to increase largely with the number of beads M because 
the phase of the estimator of 5'(g, tc) is proportional to cos <5 (see Appendix A), and it is not 
possible to give precise estimates in the limit of small Em- 

With the estimation of the accuracy of the action (for HP, we get Em = 2.5), one can easily 
determine the number of complex-time beads required in the calculation of S{q, tc) at any tc- 
M is the minimum integer for which the condition \tc\/M < e*^ is satished. Accordingly, the 
whole range of times = |A| is divided in different regions where S{q, tc) is estimated with 
a different number of beads. In practice, M = 1 for G [0, M = 2 for 2^^], 

and so on. The results obtained with this splitting are reported in Fig. [2] for the HP and 
6 = 7i/9. In the hgure, the vertical lines separate the different intervals [(M — l)^^, 
where S{q,tc) is calculated with the same number of beads M. The trends observed in 
this particular case are quite general. The results obtained are statistically reliable up to 
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FIG. 2. (Color online) Real and imaginary parts of S{q,tc) for HP, with 6 = vr/9, as a function 
of tm- The line is the exact expression (fTBIl and the points correspond to our QMC results. The 
vertical lines separate the results obtained with different number of beads M. Where not shown, 
error bars are smaller than the symbol size. 

a maximum time tm which decreases when the phase 5 is reduced. This feature directly 
implies that the maximum number of beads producing sound results is also reduced when 
approaching the real axis. In general, the number of beads is small but the high accuracy 
of the action used in the calculation makes the total covered time be quite large. In the 
case shown in Fig. [21 one can see that our QMC estimation is satisfactory up to M^ax = 3, 
with a total time tm = M^ax^m — '^•5- The results with M = 4 are spread around the 
exact function but with too large error bars to be used in the subsequent transform to the 
dynamic structure function S{q,u)). 

In Fig. El we show QMC results of the complex function S{q, tc) for the HP and different 
values of the phase 6, in comparison with the exact function flTB]) . When approaching the 
real axis, i.e. when 6 decreases, both the real and imaginary parts show an increase of 
their oscillatory behavior (notice that for HP, the exact S{q,t) for real time is periodic), 
but the maximum reachable value tm decreases. Therefore, there is a compromise between 
lowering 6 as much as possible and reaching times as large as possible. Our results show 
that the optimal phase for a posterior transform to the frequency domain is within the range 
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FIG. 3. (Color online) Real and imaginary parts of S{q,tc) for HP as a function of tm- The 
upper and lower lines in each panel correspond to the real and imaginary parts of the exact result 
(Eq. [T6l) . respectively. The symbols correspond to our QMC results (squares: real part; triangles: 
imaginary part). Each panel corresponds to the calculation of S{q,tc) for different values of the 
phase 5 of the complex time. Error bars are smaller than the symbol size. 


[7r/18,7r/9]. 

Proceeding in a similar way we have applied our method to the study of the correlation 
function for a particle in an AP. The results for the real and imaginary parts of S{q,tc) are 
shown in Fig. IHfor different values of the phase 6 ranging from 7r/2 (imaginary time) down 
to 7r/36. Our Monte Carlo results are compared with exact ones obtained by numerical 
integration over the eigenstates of the Hamiltonian (differently to the HP case, an analytical 
form for the S{q,tc) of the AP is not known). The QMC estimation of the complex-time 
correlation functions shows similar accuracy to the one achieved for the HP case. Similarly 
to HP, we recover for the AP the exact results up to a maximum value of the modulus 
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FIG. 4. (Color online) Real and imaginary parts of S{q,tc) for AP as a function of tm- The 
upper and lower lines in each panel correspond to the real and imaginary parts of the exact result, 
respectively. The symbols correspond to our QMC results (squares: real part; triangles: imaginary 
part). Each panel corresponds to the calculation of S{q,tc) for different values of the phase 6 of 
the complex time. Error bars are smaller than the symbol size. 

of the complex time tm- Beyond this value, which decreases with 5, the statistical errors 
grow signihcantly, making any estimation of S(g, tc) not reliable. Again, a good compromise 
between statistical fluctuations and approaching the real axis as close as possible locates the 
optimal values of the phase in the same range than in the HP case, 6 G [tt/IS, vr/9]. 

Once we have found the working window, the next step is to make the inversion from 
complex-time to energies. Our goal is to calculate the dynamic response S{q, oj) and compare 
the results with the exact function for both the HP and AP. To this end, we have applied 
the inversion technique described in the previous Section. A preliminary point is to know 
up to which extent the inversion procedure can influence the results in the energy domain. 
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FIG. 5. (Color online) Dynamic structure function S{q,uj) for the HP at q = 1.5. Diamonds 
correspond to the exact values and and circles and squares with errorbars to the results derived 
from the QMC results for S{q,tc)- The circles are obtained using the method described in Sec. 
II and the squares using a standard simulated annealing schedule. Left panel: imaginary time 
{5 = 7r/2). Right panel: complex time {6 = 7r/9). 


In the case of purely imaginary-time data, several inversion methods have been used,— 
the majority of them being of stochastic nature. This inverse Laplace transform is normally 
mapped to a multidimensional optimization problem. The ill-posed nature of this inversion 
can lead to results that can depend on the method employed. 

In Fig. [5l we compare results obtained for S{q,u}) in the HP problem using the inversion 
method discussed in the previous Section and a standard simulated annealing algorithm. In 
the hgure, the exact result^ 

OO 

S(q, u,) = (I’’) 

n=0 

is also plot with vertical lines. This comparison is made for two cases: imaginary-time data 
{6 = T^/^) and complex-time results with 6 = 7r/9. As it has been commented, the inversion 
from imaginary time to the frequency domain is an ill-posed problem and thus the results can 
show differences depending on the selected method. This is shown in Fig. [S](left panel): the 
inversion obtained from the stochastic simulated annealing method and the one discussed in 
Sec. II produce slightly different predictions for the higher transition lines, while they both 
agree on the hrst and second peaks, although the latter has a total strength that is ~ 15% 
off from the exact value in both cases. None of the high transition lines is well reproduced 
by any of the two models. In the same hgure (right panel), we compare the results from 
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FIG. 6. (Color online) Dynamic structure function S{q,uj) for the HP at q = 1.5. Diamonds 
correspond to the exact values (1171) and circles with errorbars to the results derived from the QMC 
results for S{q,tc). The inversion uses complex-time data calculated at the phases 5 reported in 
each panel. 


both inversion methods for 5 = vr/Q. In this case, the inversion works on complex-time 
data which shows a richer strnctnre. This signihcantly rednces the ill-posed character of the 
inversion and thus the results obtained with both methods look much more similar than in 
the 5 = 7r/2 case. Our results show that the three main peaks are well reproduced and the 
fourth one is approximated, slightly better using the non stochastic method which has been 
computed averaging over a larger data set. 

The results that we have obtained for the HP dynamic response at g = 1.5 are reported in 
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Fig. O The different panels contain reconstructions from imaginary-time data (5 = vr/2) and 


complex-time correlation factors estimated at decreasing values of the phase, down to 6 


7r/18. At 5 = 7r/2 we recover the first peak (energy and strength) and approximate the second 
one. In other words, only the lowest-energy mode is accurately reproduced. It is worth 
noticing that this is the overall trend observed in transformations from purely imaginary- 
time data. By progressively introducing a real component in the correlation factor, i.e., by 
decreasing the phase 6, the quality of the dynamic response improves significantly. As one 
can see, for 6 = n/d one gets the first four modes with their respective strengths in nice 
agreement with the exact values. By reducing even more the phase down to 5 = tt/IS we 
are able to reduce the variance of the data but no additional (higher) energies are resolved. 
Notice, however, that the strength of the peaks beyond the first four ones is much smaller. 

The same analysis has been carried out for the AP. Our results of the dynamic structure 
functions are contained in Fig. [71 With imaginary-time data, we are able to reproduce only 
the first peak. By decreasing the phase 6 the dynamic response improves progressively. At 
6 = IT/9, the three main modes and their respective strengths are in close agreement with the 
exact results. For the smallest value 6 = vr/18, we can even resolve the fourth mode whose 
strength is already quite small. Again, the gain of working with complex-time correlation 
factors becomes evident. 

When the momentum q increases, the number of modes contributing to S{q,u) also 
increases, shifting the strength to higher energies. When q is large enough, the dynamic re¬ 
sponse is centered around the free atom recoil energy ujr = h^q^/ (2m).— We have calculated 
the dynamic response for the HP and AP at g = 5. Our results are reported in Figs. [8] and 
M for HP and AP, respectively. The theoretical response shows in both cases, but somehow 
more clearly in the HP one, a distribution of modes nearly symmetric around the recoil 
energy. The results obtained for the HP are reported in Fig. [S] where we compare two cases, 
(5 = 7r/2 and S = tt/Q. Our results are shown with a continuous curve since our resolution 
does not allow for a clear separation of the individual excitation energies. Nevertheless, in 
the case of using complex time {6 = vr/9) the curve precisely reproduces the envelope of the 
exact spectrum plotted as vertical lines of strength hi given by 
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FIG. 7. (Color online) Dynamic structure function S{q,uj) for the AP at q = 1.5. Diamonds 
correspond to the exact values and circles with errorbars to the results derived from the QMC 
results for S{q,tc). The inversion uses complex-time data calculated at the phases 5 reported in 
each panel. 


located at the exact frequency modes cuj, with Aui = (cuj+i — a;j_i)/2. Using just imaginary 
time produces results which are signihcantly worse. Similar conclusions are drawn from the 
results for the AP reported in Fig. |9l The results at 5 = 7r/2 are able only to localize the 
signal of S{q,u)) around ur, but they cannot reproduce the shape of the spectral function. 
On the other hand, our results at 6 = 7r/9 match almost perfectly the exact dynamic 
response. 
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FIG. 8. (Color online) Dynamic structure function S{q,u) for the HP at q = 5. Diamonds 
correspond to the exact values dm) and the curve corresponds to the results derived from our 
QMC results. Left panel: imaginary time (d = 7r/2). Right panel: complex time (d = vr/Q). 
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FIG. 9. (Color online) Dynamic structure function S{q,uj) for the AP at q = 5. Diamonds 
correspond to the exact values and the curve to the results derived from the QMC results for 
S{q,tc)- Left panel: imaginary time {5 = 7r/2). Right panel: complex time (<5 = 7r/9). 


IV. CONCLUSIONS 


The goal of this work is to propose a new QMC strategy aimed at the study of the dynamic 
response of quantum systems at zero temperature. In quantum Monte Carlo methods, the 
evolution of configurations is carried out in purely imaginary time, both at zero and finite 
temperatures, in an attempt to describe the main properties of quantum systems with high 
accuracy. Unfortunately, dynamics in real time is not accessible and the usual approach to 
get information on the dynamic response has been to reconstruct it from purely imaginary- 
time correlation factors. However, the ill-posed character of the inverse Laplace transform 
of noisy data makes this procedure quite uncertain and with multiple solutions. 
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Our work is an attempt of progressing in a different way, that is, to reduce the ill-posed 
nature of the process by inverting data containing more information than the smooth signal 
observed in imaginary-time. Working in the zero-temperature limit, where quantumness is 
unavoidable, we have devised a strategy based on the PIGS method to sample complex-time 
correlation factors. Our method consists in the sampling of paths connecting conhgurations 
distributed according to the ground-state wave function and, in particular, the calculation 
of the correlation function in complex time over the sampled paths. The use of high-order 
actions for the propagation in complex time has proven to be crucial to get reliable data 
within a time window which naturally shrinks when the real axis is approached. Optimizing 
the phase 6 of the complex time, we have shown that, in the two model problems studied, 
we are able to improve signihcantly the calculated dynamic structure factor S{q,u)). Both 
at low and high q the description of the dynamics is signihcantly improved in comparison 
with the usual imaginary time approach. Nevertheless, additional effort is needed to conhrm 
the usefulness of the proposed method to problems in two and three dimensions and with 
more particles. Work is in progress in our group to extend this formalism to many-particle 
systems. 
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Appendix A 

In this Appendix, we report the explicit expressions for the estimator O^(xo,..., xm) ap¬ 
pearing in Eq. |H]using different actions and having chosen to sample the paths {xi, X 2 , ■ ■ ■, xm-i} 
with Ppath(a^O) Xi,..., Xm) dehned in Eq. 01 In general, Oa is a complex number that can be 
rewritten in the form 


Oa(xo, ..., Xm) 


Ylk=l^i^k,Xk-l,£c) 

Ppath(^0: ^1; • • • ; ^m) 


n 


^k—li ^c) 
^free(^/c 5 ^/c—17 


exp{C) exp{iA) , (Al) 


21 




with Ec = The terms C and A are respectively the logarithm of the modulus and the 

phase of the complex number Oa, and their formula depends on the approximation scheme 
chosen for the complex-time propagator. 

In the primitive action (PA) approximation, introducing the propagator Gpa in Eq. lAll 
we get 


M r 

CpA = 

k=l 


{xk - Xk-iY f sm6 1\ V{xk)+ V{xk+i) . „ 

-- -Sm -VT- sm6 

V ^m. 'Tq j 2/Z 


4A 


and 


M r 


Apa - Y 


k=l 


'{xk-Xk-if ^ V{xk)+ V{xk+i) ^ 

- cos 0 — Err, -- COS 0 


4:XSr 


2h 


(A2) 


(A3) 


In Chin’s approximation (CA) the propagator is given by 


Gca = JJ exp (i 

3=0 


.{Xk,j+i - Xkjf\ f y{xk,j) + V{xkj+i) 

I t 


AXtjEc 


2h 


VjEc 


X exp i 


.UoW{Xk,j) + W{Xk,j+l) 3 


C I ’ 


(A4) 


3 2h 

with a generalized potential W (rb due to the double commutator [E, [iC, E]], and parameters 
tj, Vj, and uq reported in Ref. l37|. Introducing this propagator in Eq. lAll we can End the 
functions Cca and Aca, 


4 r 

Cca = Y 
i=i 


{xkj+i -Xk,jy\ _ Y\ + 


4 At,- 


V 




^mVj 


V (Xkj+i) + V (Xkj 


2h 


A- ) sin 5 


UoW{xk,j+i) + W{xk,j) 


2h 


sin(3(5) 


(AS) 


and 


Aca - Y 
t=i 


AXtjSjYi 


cos 5+ 


V{xkj+i) + V{xkj] 

2h 


-Err,V 


COS 5 + 




^ (A6) 

3 2n I ^ > V ; 

Unfortunately, for 5 < 7r/3, the term with in the expression of Cca HAS!) is positive, and 
then exp(CcA) can become exceedingly large and spoil the calculation. 
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In order to circumvent this problem, we have worked with the sixth-order expansioii^ 


gEc-ff ^ 


64 


45 




.lg£cV'/4g£,it/2g£,V/2g£ek-/2g£ch/4 J_g£,\//2g£,h'g£,V/2 

9 45 


(A7) 


which is built without double-commutator terms. This expansion corresponds to a linear 
combination of expansions approximated with PA over the same time Sc but with different 
time steps (precisely, ec/4 in the hrst term, 6^/2 in the second term and Sc in the third 
term). Therefore, the complete formula for the exponent Cza and for the phase Aza in the 
Zillich approximation are easily obtained from Cpa and Apa (Eqs. IA2I and IA3p calculated 
for different values of Sc- 


Appendix B 

We discuss in this Appendix the method that we have followed to hnd the optimal 
regularization parameter (see Sec. II.B). Given the spectral function obtained 

inverting a series of QMC data for the complex-time correlation function C'QMc(tc) with 
a certain regularization parameter a, we calculate the complex-time correlation function 
GiNv(A,n) obtained from the integral transform of S'inv(i^,o): 

C'iNv(A,a) = J due Si-!<iy{u, a) . (Bl) 

Then we calculate the residual between C'qmc(A) and Ginv(A,o) as a function of the 
regularization parameter a. When a is large, the regularization procedure modihes the 
inversion process up to the point that GiNv(A)a) starts to differ from the previous Monte 
Carlo data C'qmc(A )5 thus showing an increase in For very small a, the noise in the 
Monte Carlo data is largely amplihed and the inversion procedure itself starts to produce 
meaningless results, giving rise once again to the increase in A plot of the total residual 
versus the regularization parameter a shows a minimum, as shown in Fig. [10] for the case 
of the AP data at g = 1.5 and 5 = 7r/4. 

In the best scenario, with high quality Monte Carlo data, an optimal regularization 
parameter may allow avoiding both problems. In any case, the full inspection of the inversion 
landscape for several values of the regularization parameter is a quick calculation. 
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FIG. 10. Residual between Cqmc(^c) and CinvC^co) as a function of the regularization param¬ 
eter a. The data corresponds to the calculation of the density correlation function S{q,tc) (Eq. 
fT^ in complex time tc = for the AP at q = 1.5 and 6 = vr/d. 
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